## Replication code for "The Conditional Effects of Public Opinion and Judicial Selection in New Judicial Federalism Cases"
## Daniel J. Mallinson and Michael Zimmerman
## 2017-03-28
## 2017-08-29
## 2017-11-28
## 2018-01-23
## 2018-02-14
## 2022-05-10

############ Reproduce and Extend Results from Zcshirnt 2009 #################
rm(list = ls()) #clears the workspace for R

#Load required libraries
library(foreign)
library(lmtest)
library(interplot)
library(grid)
library(gridExtra)
library(effects)
library(reshape)

zschirnt.data <- read.csv("zschirnt_replication_dataset.csv") #Load Zschirnt's replication data

alt.data <- read.csv("zschirnt_replication_dataset_alt_selection.csv") #Load replication data with CA, MD, and UT marked as Missouri Plan states

## Add Kastellec Dynamic MRP Estimates of State Gay Marriage Opinion
kastellec <- read.csv("state_preds_means_quantiles_from_STAN.csv")
head(kastellec)

merged.data <- merge(alt.data, kastellec, by=c("state", "year"), all=FALSE)

####################################################################################
###################################  Figure 1 ######################################
####################################################################################

## Create time series plots of state opinion (gay marriage)
kastellec.opinion.rotate <- kastellec
kastellec.opinion.rotate$mrp.trend.lower <- kastellec.opinion.rotate$mrp.trend.upper <- kastellec.opinion.rotate$mrp.trend.median <- NULL
kastellec.md <- melt(kastellec.opinion.rotate, id=c("state", "year"))
kastellec.plot.data <- cast(kastellec.md, year~variable + state)
names(kastellec.plot.data) <- c("year", levels(kastellec.md$state))
kastellec.plot.data$national_avg <- rowMeans(kastellec.plot.data[,2:51])

png("figure1.png", height=4, width=8, units="in", res=800)
plot.new()
par(mar=c(2,6,1,1))
plot.window(xlim=c(1993,2015), ylim=c(10,80))
title(main="", ylab="MRP Estimated State Support \n for Gay Marriage", xlab="", font.lab=1)
for(i in 2:51){
  lines(ts(kastellec.plot.data[,i], start=1993, frequency=1), lwd=1, col="gray")
}
lines(ts(kastellec.plot.data$national_avg, start=1993, frequency=1), lwd=2, col="black")
text(1995, 46, "Massachusetts", pos=3, font=1)
text(1995, 7, "Mississippi", pos=3, font=1)
text(2012, 56, "National Average", pos=3, font=1)
axis(1, at=seq(1993, 2015, 2), labels=seq(1993, 2015, 2), font=1)
axis(2, at=seq(10,80,10), labels=c("10%", "20%", "30%", "40%", "50%", "60%", "70%", "80%"), font=1, las=2)
#minor.tick(4, 2, tick.ratio=.5)
rug(x=seq(15,75,10), ticksize=-0.01, side=2)
rug(x=seq(1994,2014,2), ticksize=-0.01, side=1)
dev.off()



####################################################################################
###############################  Table 1  ##########################################
####################################################################################

## First, let's replicate Model 1 of Zschirnt 2017
model.1 <- glm(vote ~ conservative_judge + male + age + race + unanimous + marriage + ig_sponsorship + missouri_plan + elections + term_length + docket_control + amendment_rate + conservative_state + discrimination_law + anti_amendment + opinion, data=zschirnt.data, family=binomial(link="logit"))
summary(model.1) #replicates perfectly

## Now, with CA, MD, and UT marked as Missouri Plan states, not appointment only
model.1.alt <- glm(vote ~ conservative_judge + male + age + race + unanimous + marriage + ig_sponsorship + missouri_plan + elections + term_length + docket_control + amendment_rate + conservative_state + discrimination_law + anti_amendment + opinion, data=alt.data, family=binomial(link="logit"))
summary(model.1.alt)

## Add interactions for selection
model.1.int.all <- glm(vote ~ conservative_judge + male + age + unanimous + marriage + ig_sponsorship + missouri_plan + partisan + nonpartisan + reappointment + amendment_rate + conservative_state + anti_amendment + opinion + opinion*partisan + opinion*nonpartisan + opinion*missouri_plan + opinion*reappointment, data=alt.data, family=binomial(link="logit"))
summary(model.1.int.all)
lrtest(model.1.alt, model.1.int.all)

## Re-run Zschirnt analysis with MRP opinion estimates

model.mrp.all <- glm(vote ~ partisan + nonpartisan + missouri_plan + reappointment + mrp.mean + conservative_judge + male + age + unanimous + marriage + ig_sponsorship + amendment_rate + conservative_state + anti_amendment + mrp.mean*partisan + mrp.mean*nonpartisan + mrp.mean*missouri_plan + mrp.mean*reappointment, data=merged.data, family=binomial(link="logit"))
summary(model.mrp.all)

summary(model.1.alt, model.mrp.all)

int.a <- effect(term="partisan*mrp.mean", mod=model.mrp.all, xlevels=list(partisan=c(0,1)))
int.a$variables$partisan$is.factor <- TRUE
#levels(interaction$x$mrp.mean) <- c(0, 10, 20, 30, 40, 50, 60)
a <- plot(int.a, main="", xlab="Public Support for Gay Marriage", ylab="Predicted Pro-Gay Marriage Vote", factor.names=FALSE) 
a$condlevels$partisan <- c("Other Selection Forms", "Partisan Election")

int.b <- effect(term="nonpartisan*mrp.mean", mod=model.mrp.all, xlevels=list(nonpartisan=c(0,1)))
int.b$variables$nonpartisan$is.factor <- TRUE
#levels(interaction$x$mrp.mean) <- c(0, 10, 20, 30, 40, 50, 60)
b <- plot(int.b, main="", xlab="Public Support for Gay Marriage", ylab="Predicted Pro-Gay Marriage Vote", factor.names=FALSE) 
b$condlevels$nonpartisan <- c("Other Selection Forms", "Nonpartisan Election")

int.c <- effect(term="missouri_plan*mrp.mean", mod=model.mrp.all, xlevels=list(missouri_plan=c(0,1)))
int.c$variables$missouri_plan$is.factor <- TRUE
#levels(interaction$x$mrp.mean) <- c(0, 10, 20, 30, 40, 50, 60)
c <- plot(int.c, main="", xlab="Public Support for Gay Marriage", ylab="Predicted Pro-Gay Marriage Vote", factor.names=FALSE) 
c$condlevels$missouri_plan <- c("Other Selection Forms", "Retention Election")

int.d <- effect(term="reappointment*mrp.mean", mod=model.mrp.all, xlevels=list(reappointment=c(0,1)))
int.d$variables$reappointment$is.factor <- TRUE
#levels(interaction$x$mrp.mean) <- c(0, 10, 20, 30, 40, 50, 60)
d <- plot(int.d, main="", xlab="Public Support for Gay Marriage", ylab="Predicted Pro-Gay Marriage Vote", factor.names=FALSE) 
d$condlevels$reappointment <- c("Other Selection Forms", "Reappointment")

lay <- rbind(c(1,1,2,2),
             c(3,3,4,4))
grid.arrange(a, b, c, d, layout_matrix=lay)

####################################################################################
###############################  Table S1  #########################################
####################################################################################

model.1.int.part <- glm(vote ~ conservative_judge + male + age + unanimous + marriage + ig_sponsorship + partisan + amendment_rate + conservative_state + anti_amendment + opinion + opinion*partisan, data=alt.data, family=binomial(link="logit"))
summary(model.1.int.part)

model.1.int.nonpart <- glm(vote ~ conservative_judge + male + age + unanimous + marriage + ig_sponsorship + nonpartisan + amendment_rate + conservative_state + anti_amendment + opinion + opinion*nonpartisan, data=alt.data, family=binomial(link="logit"))
summary(model.1.int.nonpart)

model.1.int.reten <- glm(vote ~ conservative_judge + male + age + unanimous + marriage + ig_sponsorship + missouri_plan + amendment_rate + conservative_state + anti_amendment + opinion + opinion*missouri_plan, data=alt.data, family=binomial(link="logit"))
summary(model.1.int.reten)

model.1.int.reappoint <- glm(vote ~ conservative_judge + male + age + unanimous + marriage + ig_sponsorship + reappointment + amendment_rate + conservative_state + anti_amendment + opinion + opinion*reappointment, data=alt.data, family=binomial(link="logit"))
summary(model.1.int.reappoint)


int.a <- effect(term="partisan*opinion", mod=model.1.int.part, xlevels=list(partisan=c(0,1)))
int.a$variables$partisan$is.factor <- TRUE
#levels(interaction$x$opinion) <- c(0, 10, 20, 30, 40, 50, 60)
a <- plot(int.a, main="", xlab="Public Support for Gay Marriage", ylab="Predicted Pro-Gay Marriage Vote", factor.names=FALSE) 
a$condlevels$partisan <- c("Other Selection Forms", "Partisan Election")

int.b <- effect(term="nonpartisan*opinion", mod=model.1.int.nonpart, xlevels=list(nonpartisan=c(0,1)))
int.b$variables$nonpartisan$is.factor <- TRUE
#levels(interaction$x$opinion) <- c(0, 10, 20, 30, 40, 50, 60)
b <- plot(int.b, main="", xlab="Public Support for Gay Marriage", ylab="Predicted Pro-Gay Marriage Vote", factor.names=FALSE) 
b$condlevels$nonpartisan <- c("Other Selection Forms", "Nonpartisan Election")

int.c <- effect(term="missouri_plan*opinion", mod=model.1.int.reten, xlevels=list(missouri_plan=c(0,1)))
int.c$variables$missouri_plan$is.factor <- TRUE
#levels(interaction$x$opinion) <- c(0, 10, 20, 30, 40, 50, 60)
c <- plot(int.c, main="", xlab="Public Support for Gay Marriage", ylab="Predicted Pro-Gay Marriage Vote", factor.names=FALSE) 
c$condlevels$missouri_plan <- c("Other Selection Forms", "Retention Election")

int.d <- effect(term="reappointment*opinion", mod=model.1.int.reappoint, xlevels=list(reappointment=c(0,1)))
int.d$variables$reappointment$is.factor <- TRUE
#levels(interaction$x$opinion) <- c(0, 10, 20, 30, 40, 50, 60)
d <- plot(int.d, main="", xlab="Public Support for Gay Marriage", ylab="Predicted Pro-Gay Marriage Vote", factor.names=FALSE) 
d$condlevels$reappointment <- c("Other Selection Forms", "Reappointment")

lay <- rbind(c(1,1,2,2),
             c(3,3,4,4))
grid.arrange(a, b, c, d, layout_matrix=lay)

####################################################################################
###############################  Table S2  #########################################
####################################################################################

model.mrp.part <- glm(vote ~ partisan + mrp.mean + conservative_judge + male + age + unanimous + marriage + ig_sponsorship + amendment_rate + conservative_state + anti_amendment + mrp.mean*partisan, data=merged.data, family=binomial(link="logit"))
summary(model.mrp.part)

model.mrp.nonpart <- glm(vote ~ nonpartisan + mrp.mean + conservative_judge + male + age + unanimous + marriage + ig_sponsorship + amendment_rate + conservative_state + anti_amendment + mrp.mean*nonpartisan, data=merged.data, family=binomial(link="logit"))
summary(model.mrp.nonpart)

model.mrp.reten <- glm(vote ~ missouri_plan + mrp.mean + conservative_judge + male + age + unanimous + marriage + ig_sponsorship + amendment_rate + conservative_state + anti_amendment + mrp.mean*missouri_plan, data=merged.data, family=binomial(link="logit"))
summary(model.mrp.reten)

model.mrp.reappoint <- glm(vote ~ reappointment + mrp.mean + conservative_judge + male + age + unanimous + marriage + ig_sponsorship + amendment_rate + conservative_state + anti_amendment + mrp.mean*reappointment, data=merged.data, family=binomial(link="logit"))
summary(model.mrp.reappoint)

int.a <- effect(term="partisan*mrp.mean", mod=model.mrp.part, xlevels=list(partisan=c(0,1)))
int.a$variables$partisan$is.factor <- TRUE
#levels(interaction$x$mrp.mean) <- c(0, 10, 20, 30, 40, 50, 60)
a <- plot(int.a, main="", xlab="Public Support for Gay Marriage", ylab="Predicted Pro-Gay Marriage Vote", factor.names=FALSE) 
a$condlevels$partisan <- c("Other Selection Forms", "Partisan Election")

int.b <- effect(term="nonpartisan*mrp.mean", mod=model.mrp.nonpart, xlevels=list(nonpartisan=c(0,1)))
int.b$variables$nonpartisan$is.factor <- TRUE
#levels(interaction$x$mrp.mean) <- c(0, 10, 20, 30, 40, 50, 60)
b <- plot(int.b, main="", xlab="Public Support for Gay Marriage", ylab="Predicted Pro-Gay Marriage Vote", factor.names=FALSE) 
b$condlevels$nonpartisan <- c("Other Selection Forms", "Nonpartisan Election")

int.c <- effect(term="missouri_plan*mrp.mean", mod=model.mrp.reten, xlevels=list(missouri_plan=c(0,1)))
int.c$variables$missouri_plan$is.factor <- TRUE
#levels(interaction$x$mrp.mean) <- c(0, 10, 20, 30, 40, 50, 60)
c <- plot(int.c, main="", xlab="Public Support for Gay Marriage", ylab="Predicted Pro-Gay Marriage Vote", factor.names=FALSE) 
c$condlevels$missouri_plan <- c("Other Selection Forms", "Retention Election")

int.d <- effect(term="reappointment*mrp.mean", mod=model.mrp.reappoint, xlevels=list(reappointment=c(0,1)))
int.d$variables$reappointment$is.factor <- TRUE
#levels(interaction$x$mrp.mean) <- c(0, 10, 20, 30, 40, 50, 60)
d <- plot(int.d, main="", xlab="Public Support for Gay Marriage", ylab="Predicted Pro-Gay Marriage Vote", factor.names=FALSE) 
d$condlevels$reappointment <- c("Other Selection Forms", "Reappointment")

lay <- rbind(c(1,1,2,2),
             c(3,3,4,4))
grid.arrange(a, b, c, d, layout_matrix=lay)

####################################################################################
###################################  Figure 3 ######################################
####################################################################################

png("figure3.png", height=8, width=7, units="in", res=600)

a <- interplot(m=model.1.int.part, var1="opinion", var2="partisan", ci=.95, point=T) + 
  xlab("(a) Partisan Elections") +
  ylab("Marginal Effect of Pro-Gay Marriage Opinion") + 
  #ggtitle("National Trend (Zschirnt 2016)" )+
  #theme(plot.title=element_text(face="bold", hjust=0.5)) +
  theme_bw() + 
  geom_hline(yintercept=0, linetype= "dashed")

b <- interplot(m=model.1.int.nonpart, var1="opinion", var2="nonpartisan", ci=.95, point=T) + 
xlab("(b) Nonpartisan Elections") +
ylab("Marginal Effect of Pro-Gay Marriage Opinion") + 
#ggtitle("National Trend (Zschirnt 2016)" )+
#theme(plot.title=element_text(face="bold", hjust=0.5)) +
theme_bw() + 
geom_hline(yintercept=0, linetype= "dashed")

c <- interplot(m=model.1.int.reten, var1="opinion", var2="missouri_plan", ci=.95, point=T) + 
  xlab("(c) Missouri Plan Retention") +
  ylab("Marginal Effect of Pro-Gay Marriage Opinion") + 
  #ggtitle("National Trend (Zschirnt 2016)" )+
  #theme(plot.title=element_text(face="bold", hjust=0.5)) +
  theme_bw() + 
  geom_hline(yintercept=0, linetype= "dashed")

d <- interplot(m=model.1.int.reappoint, var1="opinion", var2="reappointment", ci=.95, point=T) + 
  xlab("(d) Reappointment") +
  ylab("Marginal Effect of Pro-Gay Marriage Opinion") + 
  #ggtitle("National Trend (Zschirnt 2016)" )+
  #theme(plot.title=element_text(face="bold", hjust=0.5)) +
  theme_bw() + 
  geom_hline(yintercept=0, linetype= "dashed")

lay <- rbind(c(1,1,2,2),
             c(3,3,4,4))
grid.arrange(a, b, c, d, layout_matrix=lay)
dev.off()


####################################################################################
###################################  Figure 4 ######################################
####################################################################################
png("figure4.png", height=8, width=7, units="in", res=600)

a <- interplot(m=model.mrp.part, var1="mrp.mean", var2="partisan", ci=.95, point=T) + 
  xlab("(a) Partisan Elections") +
  ylab("Marginal Effect of Pro-Gay Marriage Opinion") + 
  #ggtitle("National Trend (Zschirnt 2016)" )+
  #theme(plot.title=element_text(face="bold", hjust=0.5)) +
  theme_bw() + 
  geom_hline(yintercept=0, linetype= "dashed")

b <- interplot(m=model.mrp.nonpart, var1="mrp.mean", var2="nonpartisan", ci=.95, point=T) + 
  xlab("(b) Nonpartisan Elections") +
  ylab("Marginal Effect of Pro-Gay Marriage Opinion") + 
  #ggtitle("National Trend (Zschirnt 2016)" )+
  #theme(plot.title=element_text(face="bold", hjust=0.5)) +
  theme_bw() + 
  geom_hline(yintercept=0, linetype= "dashed")

c <- interplot(m=model.mrp.reten, var1="mrp.mean", var2="missouri_plan", ci=.95, point=T) + 
  xlab("(c) Missouri Plan Retention") +
  ylab("Marginal Effect of Pro-Gay Marriage Opinion") + 
  #ggtitle("National Trend (Zschirnt 2016)" )+
  #theme(plot.title=element_text(face="bold", hjust=0.5)) +
  theme_bw() + 
  geom_hline(yintercept=0, linetype= "dashed")

d <- interplot(m=model.mrp.reappoint, var1="mrp.mean", var2="reappointment", ci=.95, point=T) + 
  xlab("(d) Reappointment") +
  ylab("Marginal Effect of Pro-Gay Marriage Opinion") + 
  #ggtitle("National Trend (Zschirnt 2016)" )+
  #theme(plot.title=element_text(face="bold", hjust=0.5)) +
  theme_bw() + 
  geom_hline(yintercept=0, linetype= "dashed")

lay <- rbind(c(1,1,2,2),
             c(3,3,4,4))
grid.arrange(a, b, c, d, layout_matrix=lay)
dev.off()

#################################################################################
########################## Abortion Case Analysis ###############################
#################################################################################

library(foreign)
abortion.data <- read.csv("abortion_cases.csv")
abortion.opinion <- read.csv("abortion_state_preds_means_quantiles_from_STAN.csv")

## Add Dynamic MRP Estimates of State Abortion Opinion
abortion.data$merge_year <- abortion.data$year
abortion.opinion$merge_year <- abortion.opinion$year
abortion.data$merge_year[abortion.data$merge_year<1995] <- 1995 #Need to use 1995 opinion data for 1994 decisions
abortion.data <- merge(abortion.data, abortion.opinion, by=c("state", "merge_year"), all=FALSE)

## Combine partisan and non-partisan elections
abortion.data$elections <- abortion.data$nonpartisan + abortion.data$partisan

## Calculate justice age at time of opinion
abortion.data$age <- abortion.data$year.x - abortion.data$birth

####################################################################################
###################################  Figure 2 ######################################
####################################################################################
## Create time series plots of state opinion
library(reshape)
abortion.opinion.rotate <- abortion.opinion
abortion.opinion.rotate$mrp.trend.lower <- abortion.opinion$mrp.trend.upper <- abortion.opinion$mrp.trend.median <- NULL
abort.md <- melt(abortion.opinion.rotate, id=c("state", "year"))
abortion.plot.data <- cast(abort.md, year~variable + state)
names(abortion.plot.data) <- c("year", levels(abort.md$state))
abortion.plot.data$national_avg <- rowMeans(abortion.plot.data[,2:51])

png("figure2.png", height=4, width=8, units="in", res=800)
plot.new()
par(mar=c(2,6,1,1))
plot.window(xlim=c(1995,2015), ylim=c(30,80))
title(main="", ylab="MRP Estimated State Support \n for Abortion Availability", xlab="", font.lab=1)
for(i in 2:51){
lines(ts(abortion.plot.data[,i], start=1995, frequency=1), lwd=1, col="gray")
}
lines(ts(abortion.plot.data$national_avg, start=1995, frequency=1), lwd=2, col="black")
text(1998, 74, "New Hampshire", pos=3, font=1)
text(1998, 33, "Mississippi", pos=3, font=1)
text(2012, 53, "National Average", pos=3, font=1)
axis(1, at=seq(1995, 2015, 5), labels=seq(1995, 2015, 5), font=1)
axis(2, at=seq(30,75,5), labels=seq(30,75,5), font=1, las=2)
rug(x=c(31:34, 36:39, 41:44, 46:49, 51:54, 56:59, 61:64, 66:69, 71:74), ticksize=-0.01, side=2)
rug(x=c(1996:1999, 2001:2004, 2006:2009, 2011:2014), ticksize=-0.01, side=1)
dev.off()

########################### Analyses ###################################

# Create two analysis datasets, one with duplicate votes and one without
abortion.data.nodup <- abortion.data[which(abortion.data$vote!=0.5),]

####################################################################################
###################################  Table 2  ######################################
####################################################################################

#No interactions - pared from Zschirnt (2016)
model.noint.nodup <- glm(vote ~ partisan + nonpartisan + missouri_plan + conservative_judge + male + age + unanimous + amendment_rate + conservative_state + mrp.mean, data=abortion.data.nodup, family=binomial(link="logit"))
summary(model.noint.nodup)
lrtest(model.noint.nodup)

# Further reduced model
model.reduced.nodup <- glm(vote ~ partisan + nonpartisan + missouri_plan + amendment_rate + mrp.mean, data=abortion.data.nodup, family=binomial(link="logit"))
summary(model.reduced.nodup)
lrtest(model.reduced.nodup)

## Individual models with selection method interactions

## Problems with separation when partisan elections added
model.intpart.nodup <- glm(vote ~ partisan + mrp.mean + amendment_rate + partisan*mrp.mean, data=abortion.data.nodup, family=binomial(link="logit"))

#Use penalized maximum likelihood to deal with separation
#install.packages("logistf")
#library(logistf)
#install.packages("brglm")
library(brglm)
model.intpart.nodup <- brglm(vote ~ partisan + mrp.mean + amendment_rate + partisan*mrp.mean, data=abortion.data.nodup, family=binomial(link="logit"))
summary(model.intpart.nodup)
lrtest(model.intpart.nodup)

model.intnon.nodup <- glm(vote ~ nonpartisan + mrp.mean + amendment_rate + nonpartisan*mrp.mean, data=abortion.data.nodup, family=binomial(link="logit"))
summary(model.intnon.nodup)
lrtest(model.intnon.nodup)

model.intreten.nodup <- glm(vote ~ missouri_plan + mrp.mean + amendment_rate + missouri_plan*mrp.mean, data=abortion.data.nodup, family=binomial(link="logit"))
summary(model.intreten.nodup)
lrtest(model.intreten.nodup)

#NJ is the only reappointment state, model doesn't tell much
#model.intappt.nodup <- glm(vote ~ reappointment + mrp.mean + amendment_rate + reappointment*mrp.mean, data=abortion.data.nodup, family=binomial(link="logit"))
#summary(model.intappt.nodup)
#lrtest(model.intappt.nodup)

####################################################################################
##################################  Figure 5  ######################################
####################################################################################

png("figure5.png", height=8, width=7, units="in", res=600)

int.a <- effect(term="partisan*mrp.mean", mod=model.intpart.nodup, xlevels=list(partisan=c(0,1)))
int.a$variables$partisan$is.factor <- TRUE
#levels(interaction$x$mrp.mean) <- c(0, 10, 20, 30, 40, 50, 60)
a <- plot(int.a, main="", xlab="Public Support for Abortion Access", ylab="Predicted Pro-Access Vote", factor.names=FALSE) 
a$condlevels$partisan <- c("Other Selection Forms", "Partisan Election")

int.b <- effect(term="nonpartisan*mrp.mean", mod=model.intnon.nodup, xlevels=list(nonpartisan=c(0,1)))
int.b$variables$nonpartisan$is.factor <- TRUE
#levels(interaction$x$mrp.mean) <- c(0, 10, 20, 30, 40, 50, 60)
b <- plot(int.b, main="", xlab="Public Support for Abortion Access", ylab="Predicted Pro-Access Vote", factor.names=FALSE) 
b$condlevels$nonpartisan <- c("Other Selection Forms", "Nonpartisan Election")

int.c <- effect(term="missouri_plan*mrp.mean", mod=model.intreten.nodup, xlevels=list(missouri_plan=c(0,1)))
int.c$variables$missouri_plan$is.factor <- TRUE
#levels(interaction$x$mrp.mean) <- c(0, 10, 20, 30, 40, 50, 60)
c <- plot(int.c, main="", xlab="Public Support for Abortion Access", ylab="Predicted Pro-Access Vote", factor.names=FALSE) 
c$condlevels$missouri_plan <- c("Other Selection Forms", "Retention Election")

#int.d <- effect(term="reappointment*mrp.mean", mod=model.intappt.nodup, xlevels=list(reappointment=c(0,1)))
#int.d$variables$reappointment$is.factor <- TRUE
#levels(interaction$x$mrp.mean) <- c(0, 10, 20, 30, 40, 50, 60)
#d <- plot(int.d, main="", xlab="Public Support for Abortion Access", ylab="Predicted Pro-Access Vote", factor.names=FALSE) 
#d$condlevels$reappointment <- c("Other Selection Forms", "Reappointment")
#d

lay <- rbind(c(1,1,2,2),
             c(4,3,3,4))
grid.arrange(a, b, c, layout_matrix=lay)
dev.off()

#png("abortion_partisan_interaction.png", height=5, width=6, units="in", res=600)
#interaction <- effect(term="partisan*mrp.mean", mod=model.intpart.nodup, xlevels=list(partisan=c(0,1)))
#interaction$variables$partisan$is.factor <- TRUE
#levels(interaction$x$mrp.mean) <- c(0, 10, 20, 30, 40, 50, 60)
#interact_plot <- plot(interaction, main="", xlab="Public Support for Abortion Access", ylab="Predicted Pro-Access Vote", factor.names=FALSE) 
#interact_plot$condlevels$partisan <- c("Other Selection Forms", "Partisan Election")
#interact_plot
#dev.off()

#####################################################################################
############################## Table S3 #############################################
#####################################################################################
model.noint <- glm(vote_duplicate ~ partisan + nonpartisan + missouri_plan + conservative_judge + male + age + unanimous + amendment_rate + conservative_state + mrp.mean, data=abortion.data, family=binomial(link="logit"))
summary(model.noint)
lrtest(model.noint)

model.reduced <- glm(vote_duplicate ~ partisan + nonpartisan + missouri_plan + amendment_rate + mrp.mean, data=abortion.data, family=binomial(link="logit"))
summary(model.reduced)
lrtest(model.reduced)

model.intpart <- glm(vote_duplicate ~ partisan + mrp.mean + amendment_rate + partisan*mrp.mean, data=abortion.data, family=binomial(link="logit"))
summary(model.intpart)
lrtest(model.intpart)

#library(interplot)
#png("interplot_abortion.png", height=6, width=6, units="in", res=800)
#interplot(m=model.intpart, var1="mrp.mean", var2="partisan", ci=.95, point=T) + 
#  xlab("States without (0) or with (1) Partisan Elections") +
#  ylab("Marginal Effect of Pro-Abortion Accessibility Opinion") + 
  #ggtitle("State-Level MRP" )+
  #theme(plot.title=element_text(face="bold", hjust=0.5)) +
#  theme_bw() + 
#  geom_hline(yintercept=0, linetype= "dashed")
#dev.off()

model.intnon <- glm(vote_duplicate ~ nonpartisan + mrp.mean + amendment_rate + nonpartisan*mrp.mean, data=abortion.data, family=binomial(link="logit"))
summary(model.intnon)
lrtest(model.intnon)

model.intmiss <- glm(vote_duplicate ~ missouri_plan + mrp.mean + amendment_rate + missouri_plan*mrp.mean, data=abortion.data, family=binomial(link="logit"))
summary(model.intmiss)
lrtest(model.intmiss)

